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Abstract. This paper presents algorithms for solving niultiobjective integer programming 
problems. The algorithm uses Barvinok's rational functions of the polytope that defines the 
feasible region and provides as output the entire set of nondominated solutions for the problem. 
Theoretical complexity results on the algorithm are provided in the paper. Specifically, we prove 
that encoding the entire set of nondominated solutions of the problem is polynomially doable, 
when the dimension of the decision space is fixed. In addition, we provide polynomial delay 
algorithms for enumerating this set. An implementation of the algorithm shows that it is useful 
for solving multiobjective integer linear programs. 

Introduction 

Short rational functions were used by Barvinok [2J as a tool to develop an algorithm for counting 
the number of integer points inside convex polytopes, based in the previous geometrical papers 
by Brion [6] ,Khovanskii and Puhlikov [22], and Lawrence [21]. The main idea is encoding those 
integral points in a rational function in as many variables as the dimension of the space where 
the body lives. Let P C M'^ be a given convex polyhedron, the integral points may be expressed 
in a formal sum /(P, z) = J2a with a = (ai, . . . , a^) G P fl Z'^, where 2" = z^^ ■ ■ ■ z^"^. 
Barvinok's aimed objective was representing that formal sum of monomials in the multivariate 
polynomial ring Z[zi, . . . , z^], as a "short" sum of rational functions in the same variables. 
Actually, Barvinok presented a polynomial-time algorithm when the dimension, n, is fixed, to 
compute those functions. A clear example is the polytope P = [0, A^] C M: the long expression 
of the generating function is /(P, z) = J2iLo ^'^ ■> ^"^^ ^^^y Xhi^i its representation as 

sum of rational functions is the well known formula . 

Brion proved in 1988 [6|, that for computing the short generating function of the formal sum 
associated to a polyhedron, it is enough to do it for tangent cones at each vertex of P. Barvinok 
applied this function to count the number of integral points inside a polyhedron P, that is, 
lim^^(i 1) /(P, 2;), that is not possible to compute using the original expression, but it may be 
obtained using tools from complex analysis over the rational function /. 

The above approach, apart from counting lattice points, has been used to develop some algo- 
rithms to solve, exactly, integer programming. Actually, De Loera et al [9| and Woods and 
Yoshida [35] presented different methods to solve this family of problems using Barvinok's 
rational function of the polytope defined by the constraints of the given problem. 
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The goal of this paper is to present new methods for solving multiobjective integer programming 
problems. In contrast to usual integer programming problems, in multiobjective problems there 
are at least two (and in this case the problem is called biobjective) or more objective functions 
to be optimized. 

The importance of multiobjective optimization is not only due to its theoretical implications 
but also to its many applications. Witnesses of that are the large number of real-world deci- 
sion problems that appear in the literature formulated as multiobjective programs. Examples 
of them are flowshop scheduling (see [IB]), analysis in Finance (see [13], Chapter 20), rail- 
way network infrastructure capacity (see [H]), vehicle routing problems (see [201 129] ) or trip 
organization (see [3T]). 

Multiobjective programs are formulated as optimization (we restrict ourselves without lost of 
generality to the maximization case) problems over feasible regions with at least two objective 
functions. Usually, it is not possible to maximize all the objective functions simultaneously 
since objective functions induce a partial order over the vectors in the feasible region, so a 
different notion of solution is needed. A feasible vector is said to be a nondominated (or Pareto 
optimal) solution if no other feasible vector has componentwise larger objective values. The 
evaluation through the objectives of a nondominated solution is called efficient solution. 
This paper studies multiobjective integer linear programs (MOILP). Thus, we assume that there 
are at least two objective functions involved, the constraints that define the feasible region are 
linear, and the feasible vectors are integers. 

Even if we assume that the objective functions are also linear, there are nowadays relatively 
few exact methods to solve general multiobjective integer and linear problems (see [13]). Some 
of them, as branch and bound with bound sets, which belong to the class of implicit enumera- 
tion methods, combine optimality of the returned solutions with adaptability to a wide range 
of problems (see for example [3Sl EH ESI EZ] for details). Other methods, as Dynamic Pro- 
gramming, are general methods for solving, not very efficiently, general families of optimization 
problems (see [211 C|)- A different approach, as the Two-Phase method (see [33]), looks for 
supported solutions (those that can be found as solutions of a single-objective problem over 
the same feasible region but with objective function a linear combination of the original ob- 
jectives) in a first stage and non-supported solutions are found in a second phase using the 
supported ones. The Two-phase method combines usual single-criteria methods with specific 
multiobjective techniques. 

Apart from those generic methods, there are specific algorithms for solving some combinatorial 
biobjective problems: biobjective knapsacks ([3l|), biobjective minimum spanning tree prob- 
lems ([30]) or biobjective assignment problems ([28]), as well as heuristics and metaheuristics 
algorithms that decrease the CPU time for computing the nondominated solutions for specific 
biobjective problems. 

Nowadays, new approaches for solving multiobjective problems, using tools from Algebraic 
Geometry and Computational Algebra, have been proposed in the literature aiming to provide 
new insights into the combinatorial structure of the problems. This new research line seems to 
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be prolific in a near future. An example of that is presented in [5] where a notion of partial 
Grobner basis is given that allows to build a test family (analogous to the test set concept but 
for solving multiobjective problems) to solve general multiobjective linear integer programming 
problems. 

Another witness of this trend is the recent work by Deloera et al. [10]. In this paper, the 
authors present several algorithms for multiobjective integer linear programs using generating 
functions. Nevertheless, their approach differs from ours in that their requires, in addition, 
to fix the dimension of the objective space to prove polynomiality of their algorithms, and 
their proofs are totally different. Moreover, no actual implementation of the algorithms is 
shown in that paper although it addresses an interesting shortest distance problems respect to 
a prespecified Pareto point. 

In this paper, we also use rational generating function of polytopes for solving multiobjective 
integer linear programs. 

In Section 1, the main results on Barvinok's rational functions, which we use in our approach, are 
presented. Section 2 presents the multiobjective integer problem and the notion of dominance 
in order to clarify which kind of solutions we are looking for. The two following sections 
analyze different algorithms for solving general multiobjective problems. In Section 3, fixing 
the dimension of the decision space, a polynomial time algorithm that encodes the set of 
nondominated solutions of the problem as a short sum of rational functions is detailed. Next, a 
digging algorithm that computes the entire set of nondominated solutions using the multivariate 
Laurent expansion for the Barvinok's function of the polytope defined by the constraints of 
the problem is given in Section 4. In that section, a polynomial delay algorithm for solving 
multiobjective problems is also presented (fixing only the dimension of the decision space). 
Section 5 shows the results of a computational experiment and its analysis. Here, we solve 
biobjective knapsack problems, report on the performance of the algorithms and draw some 
conclusions on their results and their implications. 

1. Barvinok's rational functions 

In this section, we recall some results on short rational functions for polytopes, that we use in 
our development. For details the interested reader is referred to [21 |3l H]. 

Let P = {x G M" : Ax < 6} be a rational polytope in M". The main idea of Barvinok's Theory 
was encoding the integer points inside a rational polytope in a "long" sum of monomials: 

f{P;z)= 

where = ■ ■ ■ z"". 

The following results, due to Barvinok, allow us to re-encode, in polynomial-time for fixed 
dimension, these integer points in a "short" sum of rational functions. 

Theorem 1.1 (Theorem 5.4 in [2])- Assume n, the dimension, is fixed. Given a rational 
polyhedron P C , the generating function f{P; z) can he computed in polynomial time in the 
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form 

mz) = J2e,- 

where I is a polynomial- size indexing set, and where e G {1, —1} and Ui,Vij G Z" for all i and 
j- 

As a corollary of this result, Barvinok gave an algorithm for counting the number of integer 
points in P. It is clear from the original expression of f{P;z) that this number is f{P; 1), 
but 1 = (1, . . . , 1) is a pole for the rational function, so, the number of integer points in the 
polyhedron is lim/(S';2;). This limit can be computed using residue calculation tools from 
elementary complex analysis. 

Another useful result due to Barvinok and Wood [4J, states that computing the short rational 
function of the intersection of two polytopes, given the respective short rational function for 
each polytope, is doable in polynomial time. 

Theorem 1.2 (Theorem 3.6 in [4J). Let Pi, P2 be polytopes in M" and P = PinP2- Let /(Pi; z) 
and f{P2',z) he their short rational functions with at most k binomials in each denominator. 
Then there exists a polynomial time algorithm that computes 

f{P■,z) = J2l^- 

fj(l-z''»0 

i=i 

with s < 2k, where the 7^ are rational numbers and Ui, Vij are nonzero integral vectors for i ^ I 
and j = 1, . . . , s. 

In the proof of the above theorem, the Hadamard product of a pair of power series is used. 
Given gi{z) = (3m and g2{z) = 7m -z™, the Hadamard product g = gi * g2 is the 

power series 

9iz) = ^ Vmz"^ where rjm = Pmlm- 



The following Lemma is instrumental to prove Theorem 1.2 



Lemma 1.1 (Lemma 3.4 in |3]). Let us fix k. Then there exists a polynomial time algorithm, 
which, given functions gi{z) and g2{z) such that 

(1) gi{z) = — — — and g2{z) 



1 _ ^aiij . . . (1 _ ^aifej ' (l_ ^aaij . . . (1 _ z'^^k) 

where Pi,aij G Z"' and such that there exists I G with {l,aij) < for all i,j, computes a 
function h{z) in the form 



(1 _ ^fe.i) . . . (1 _ 2^'.^) 
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with qi, bij G Z'^, /3j G Q and s < 2k such that h possesses the Laurent expansion in a neighbor- 
hood U of zq = (e'l, . . . , e'") and h{z) = gi{z) * g2{z). 



For proving Theorem 1.2, it is enough to assure that for given polytopes Pi,P2 C Z", their 
rational functions satisfy conditions ([T]). It is not difficult to ensure that the conditions are 
verified after some changes are done in the expressions for the short rational functions (for 
further details, the interested reader is referred to 

Actually, with this result a general theorem can be proved ensuring that for a pair of polytopes. 
Pi, P2 ^ there exists a polynomial time algorithm to compute, given the rational functions 
for Pi and P2, the short rational function of any boolean combination of Pi and P2. 
Finally, we recall that one can find, in polynomial time, rational functions for polytopes that 
are images of polytopes with known rational function. 

Lemma 1.2 (Theorem 1.7 in |2j). Let us fix n. There exists a numbers = s{n) and a polynomial 
time algorithm, which, given a rational polytope P C M" and a linear transformation T : M" 
W such that T(Z") C computes the function f{S;z) for S = T{P n Z"), S C Z^ zn the 
form 

" 1^ "7l-;2«-)---(l-^«»«) 

where ai G Q, pi, aij G Z^ and aij 7^ for all i,j. 

To finish this section, we mention the application of short rational functions to solve single- 
objective integer programming. The interested reader is referred to [HI E5] for further details. 

2. MULTIOBJECTIVE COMBINATORIAL OPTIMIZATION PROBLEMS 

In this section we present the problem to be solved as well as the new concept of solutions 
motivated by the nature of the problem. 

A multiobjective integer linear program (MOILP) can be formulated as: 

max (ci X, . . . , Cfc x) =: C x 
s.t. 

d 

(2) aij Xj <bi i = 1, . . . ,m 

i=i 

Xj G Z_|_ j = 1, . . . ,n 



with ttij, bi integers and Xi non negative. Without loss of generality, we will consider the above 
problem in its standard form, i.e., the coefficient of the k objective functions are non-negative 
and the constraints are in equation form. In addition, we will assume that the constraints 
define a polytope (bounded) in M". Therefore, from now on we deal with MOILPA,c{b)- 
It is clear that Problem ([2]) is not a standard optimization problem since the objective function is 
a /c-coordinate vector, thus inducing a partial order among its feasible solutions. Hence, solving 
the above problem requires an alternative concept of solution, namely the set of nondominated 
(or Pareto-optimal) points. 
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A vector x G M" is said to be a nondominated (or Pareto optimal) solution of MIOLPa,c if 
there is no other feasible vector y such that 

Cjy>CjX Vj = l,...,/c 

with at least one strict inequality for some j. If x is a nondominated solution, the vector 
Cx = (ci X, . . . , Cfc x) G M.'' is called efficient. Note that Xe is a subset of MJ^ {decision space) 
and Ye is a subset of M.'^ {objectives space). 

We will say that a dominated point, y, is dominated by x if q x ^ q ?/ for all i = 1, . . . , 
We denote by Xe the set of all nondominated solutions for ^ and by Ye the image under the 
objective functions of Xe, that is. Ye = {C x : x G Xe}- 

From the objective function C, we obtain a linear partial order on Z" as follows: 

X >-c y Cx ^ 

Notice that since C G Z*^^", the above relation is not complete. Hence, there may exist 
non-comparable vectors. We will use this partial order, induced by the objective functions of 
problem (|2| as the input for the multiobjective integer programming algorithm developed in 
this paper. 

Sometimes, the same efficient value is the image of several nondominated solutions. At this 
point, different problems can be tackled. We say that two nondominated solutions, x^ and x^ 
are equivalent if C x^ = C x^. Then, the solutions for MOILPA,c{b) are one of the following: 

• Complete set: A subset X C Xe such that for all y eYe there is x G X with C x = y. 

• Minimal complete set: A complete set with no equivalent solutions. 

• Maximal complete set: All equivalent solutions. 

Through this paper, we are looking for the entire set of nondominated solutions, equivalently 
the maximal complete set for MOILPa^- 



3. A SHORT RATIONAL FUNCTION EXPRESSION OF THE ENTIRE SET OF NONDOMINATED 

SOLUTIONS 

We present in this section an algorithm for solving MOILPAfi{h) using Barvinok's rational 
functions technique. 

Theorem 3.1. Let A G Z™^", 6 G Z™, C = (ci, . . . , c^) G Z'^^", J G {1, . . . , n}, and assume 
that the number of variables n is fixed. Suppose P = {x G : Ax < 6, x > 0} is a rational 
convex polytope in M". Then, we can encode, in polynomial time, the entire set of nondominated 
solutions for MOILPA,c{b) in a short sum of rational functions. 



We are denoting by ^ the binary relation "more than or equal to" and where it is assumed that at least one 
of the inequalities in the list is strict. 
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Proof. Using Barvinok's algorithm, compute the following generating function in 2n variables: 



(3) fix,y):= Yl 

(w,'u)ePcnz2" 

k 

where Pc = {{,u,v) G Z" x Z" : u,v G P,CiU — CiV > for alH = 1, . . . , /c and CiU — 



1=1 



^^Cjf > 1}. Pc is clearly a rational polytope. For fixed u G Z", the ^/-degrees, a, in the 

i=l 

monomial x"?/" of f{x,y) represent the solutions dominated by u. 

Now, for any function ip, let ni,ip, T^2,ip be the projections of ip{x, y) onto the x- and y-variables, 
respectively. Thus 1^2 j{y) encodes all dominated feasible integral vectors (because the degree 
vectors of the x- variables dominate them, by construction), and it can be computed from / (a;, y) 



in polynomial time by Lemma 1.2 



Let V{P) be the set of extreme points of the polytope P and choose an integer R > ma.x{vi : 
V G V{P),i = 1,...,?T,} (we can find such an integer R via linear programming). For this 
positive integer, R, let r{x, R) be the rational function for the polytope {u G : Ui < R}, its 
expression is: 



r(x,i?) = nf-^ + -^) 



Define f{x,y) as above, vr2j(x) the projection of / onto the second set of variables as a func- 
tion of the x-variables and F{x) the short generating function of P. They are computed in 
polynomial time by Lemma [L2] and Theorem |Ll| respectively. Compute the following difference: 

h{x) := F{x) — iT2j{x). 

This is the sum over all monomials where m G P is a non dominated solution, since we are 
deleting, from the total sum of feasible solutions, the set of dominated ones. 
This construction gives us a short rational function associated with the sum over all monomials 
with degrees being the nondominated solutions for MOILPA^cip). As a consequence, we can 
compute the number of nondominated solutions for the problem. The complexity of the entire 
construction being polynomial since we only use polynomial time operations among four short 
rational functions of polytopes (these operations are the computation of the short rational 
expressions for f{x,y), r{x,R) and 7r2j(a;)). □ 

Remark 3.1. To prove the above result one can use a different approach to compute the non- 
dominated solutions assuming that there exists a polynomially bounded (for fixed dimension ) 
feasible lower bound set, L, for MOILPA^ci^), i-^-, o set of feasible solutions such that every 
nondominated solution is either one element in L, or it dominates at least one the elements in 
L. 

First, compute the following operations with generating functions: 

H{x, y) = f{x, y) - /(x, y) * (vr2j(x) r{y, R)) 
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This is the sum over all monomials where u,v G P, u is a nondominated solution and v 
is dominated by u. In H{x,y), each nondominated solution, u, appears as many times as the 
number of feasible solutions that it dominates. 

Next, compute a feasible lower bound set (see [161 EJA L = . . . ,as}- This way the set of 
nondominated solutions is encoded using the following construction: 
Let RLB^{x, y) be the following short sum of rational functions 

RLB\x, y) = H{x, y) * (y"' r(x, R)) z = l,...,s. 

Taking into account that for each i , the element y""^ is common factor for RLB^{x,y) and it is the 

RLB^ix ij) 

unique factor where the y-variables appear, we can define ND^(x) = ^ — , i = 1, . . . , s, 

to be the sum of rational functions that encodes the nondominated solutions that dominate ai, 
2 = 1, . . . , s. Therefore, the entire set of nondominated solutions for MOILP^^cib) is encoded 

k 

in the short sum of rational functions ND(x) = ND^{x). 

i=l 



4. Digging algorithm for the set of nondominated solutions of MOILP 

Section 3 proves that encoding the entire set of efficient solutions of MOILP can be done in 
polynomial time for fixed dimension. This is a compact representation of the solution concept. 
Nevertheless, one may be interested in an explicit description of this list of points. This task 
could be performed, by expanding the short rational expression which is ensured by Theorem 



3.1 but it would require the implementation of all operations used in the proof. As far we 
know, they have never been efficiently implemented. 

An alternative algorithm for enumerating the nondominated solutions of a multiobjective inte- 
ger programming problem, which uses rational generating functions, is the digging algorithm. 
This algorithm is an extension of a heuristic proposed by Lasserre [23] for the single-objective 
case. 

Let A,C and b be as in Problem ([2]), and assume that P = {x E M" : Ax < b,x > 0} 
is a polytope. Then, by Theorem |1.1[ we can compute a rational expression for /(P; z) = 
EaePnZ" ^" the form 

f{P;z) = J2e,- 

JJ(l-;z'''0 

in polynomial time for fixed dimension, n. Each addend in the above sum will be referred to 
as fi, i G /. 

If we make the substitution Zi = Zi t^^' ■ ■ ■ t^*"' , in the monomial description we have /(P; z,ti, . . . ,tk) 
z''tl"' ■■■tl'", where ci, . . . , Cfc are the rows in C. It is clear that for enumerating the 

entire set of nondominated solutions, it would suffice to look for the set of leader terms, in the 
i-variables, in the partial order induced by C, >~c, of the multi-polynomial /(P; z,ti, . . . ,tk). 
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After the above changes we have: 

(4) /(P; z,h,..., tk) = f'^^'^ ^^ti,..., tk), 

where fi{P; z,ti, . . . ,tk) Si— ^ . Now, we can assume, wlog, that ci Vij 

is negative or zero. If it were zero, then we could assume that C2Vij is negative. Otherwise, 
we would repeat the argument until the first non zero element is found (it is assured that this 
element exists, otherwise the factor would not appear in the expression of the short rational 
function). Indeed, if the first non zero element were positive, we would make the change: 



1 - Z'"^o fp' ■ ■ ■ ff'' 1 - z-''^^ ti"''''' ■ ■ ■ tl"'''''' 
and the sign of the ti-degree would be negative. 

With these assumptions, the multivariate Laurent series expansion for each rational function, 
/i, in f{P;z,ti,...,tk) is 

d oo d 
3 = 1 A=0 3 = 1 

The following result allows us to develop a finite algorithm for solving MOILPA,c{b) using 
Barvinok's rational generating functions. 



Lemma 4.1. Obtaining the entire set of nondominated solutions for a MOILP requires only 
an explicit finite, polynomially bounded (in fixed dimension) number of terms of the long sum 
in the Laurent expansion of f{P; z,ti, . . . , tk). 

n 

Proof. Let i e /, j e {1, . . . , n} and define Pj = (A e Z" : CgUi + XrCsVir > 0,s = 1, . . . , k}, 

r=l 

Mij — max{Aj : A e Pj} and m^- = min{Aj : A e Pj}. My and rriij are well-defined because Pj, 
defined above, is non empty and bounded since, by construction, for each j e {1, . . . ,n} there 
exists s & {1, . . . ,k} such that Vij < 0. 

Then, it is enough to search for the nondominated solutions in the finite sum 

d Mij 
j=l \=mij 

Let U (resp. /) be the greatest (resp. smallest) value that appears in the non-zero absolute 
values of the entries in A, b, C. Set M = max{[/, /^^}. First, rriij > 0. Then, by applying 
Cramer's rule one can see that Mij is bounded above by ©(M^""*"^). Thus, the explicit number 

n 

of terms in the expansion of fi, namely — rriij \, is polynomial, when the dimension, n 

is fixed. □ 
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The digging algorithm looks for the leader terms in the t- variables, with respect to the partial or- 
der induced by C . At each rational function (addends in the above sum Q) multiplications are 
done in lexicographical order in their respective bounded hypercubes. If the t-degree of a specific 
multiplication is not dominated by one of the previous factors, it is kept in a list; otherwise the 
algorithm continues augmenting lexicographically the lambdas. To simplify the search at each 
addend, the following consideration can be taken into account: if t^^^^ ^ . . -t^"^^^ ^ ^ is 
dominated, then any term of the form ^^^^ ■ ■ ■ ^^"j ^ ^ being componentwise larger 

than A, is dominated as well. 

The above process is done on each rational function that appears in the representation of /. As 
an output we get a set of leader terms (for each rational function), that are the candidates to be 
nondominated solutions. Terms that appear with opposite signs will be cancelled. Removing 
terms in the list of candidates (to be nondominated solutions) implies consideration of those 
terms that were dominated by the cancelled ones. These terms are included in the current list 
of candidates and the process continues until no more terms are added. 

At the end, some dominated elements may appear in the union of the final list. Deleting them 
in a simple cleaning process gives the list that contains only the entire set of nondominated 
solutions for the multiobjective problem. 
Algorithm [T] details the pseudocode of the digging algorithm. 

Recall that M = max{?7, /~^}, where U is the greatest value that appears in the non-zero 
absolute values of the entries in A, b, C and / is the least value among these values. 
Taking into account Lemma 4A_ and the fact that Algorithm [T] never cycles, we have the following 
statement. 

Theorem 4.1. Algorithm^ computes in a finite (bounded on M) number of steps, the entire 
set of nondominated solutions of the multiobjective Problem 



It is well known that enumerating the nondominated solutions of MOILP is NP-hard and #P- 
hard ([ElIIl]). Thus, one cannot expect to have very efficient algorithms for solving the general 
problem (when the dimension is part of the input). 

In the following, we concentrate on a different concept of complexity that has been already 
used in the literature for slightly different problems. Computing maximal independent sets on 
graphs is known to be #P-hard ([17j), nevertheless there exist algorithms for obtaining these 
sets which ensure that the number of operations necessary to obtain two consecutive solutions 
of the problem is bounded by a polynomial in the problem input size (see e.g. [32] )• These 
algorithms are called polynomial delay. Formally, an algorithm is said polynomial delay if the 
delay, which is the maximum computation time between two consecutive outputs, is bounded 
by a polynomial in the input size ([H [T9]). 

In our case, a polynomial delay algorithm, in fixed dimension, for solving a multiobjective 
linear integer program means that once the first nondominated solution is computed, either in 
polynomial time a next nondominated solution is found or the termination of the algorithm is 
given as an output. 
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Algorithm 1: Digging algorithm for multiobjective problems 



input : A e Z™^", 6 G Z", C G Z'^^" 
Step 1: (Initialization) 

Compute, f{z), the short sum of rational functions encoding the set of nondominated solutions of 
MOILPA,c{b). The number of rational function is indexed by /. 

Make the substitution Zi = zi t^^' ■ ■ ■ t''^^ in f{z). Denote by fi, i E I, each one of the addends in 
/, as in Q. 

Set rriij and Mij, j = 1, . . . ,n, the lower and upper bounds computed in the proof of Lemma 4.1 



and S = Y\[iTT'ij , Mij] fl Z". Set Fj := {}, i e I, the initial set of nondominated solutions 

i=i 

encoded in /j. 

Step 2: (Nondominance test) 

repeat 

for i G / do 

for A* G 5 such that its entries are not componentwise larger than a previous A do 

n 

Compute Pi := z"^" t^^ ■ ■ ■ t^'' , being Wo := Ui + ^ A* Vij and 

i=i 

n 

Wh ■■= ciUi+ ^ A} ChVij h= 1,. . . ,k 
i=i 

if p is nondominated by elements in Fj then Fj ^ Fj U {p} 
end 
end 

Step 3: (Feasibility test) 

for s,r E I, s < r do 

I if p G F, n F;,, ej = -Eh then F, ^ F, \ {p}; Th^Vh\ {p} 
end 

until No changes in any Fj are done for all i E I ; 

Set F := [JFj. Remove from F the dominated elements. 

j 

output: The entire set of nondominated solutions for MOILPa,c{^)'- T 



Next, we present a polynomial delay algorithm, in fixed dimension, for solving multiobjective 
integer linear programming problems. This algorithm combines the theoretical construction of 



Theorem 3.1 and a digging process in the Laurent expansion of the short rational functions of 



the polytope associated with the constraints of the problem. 



The algorithm proceeds as follows. 
Let f{z) be the short rational function that encodes the nondominated solutions (by Theorem 



3.1, the complexity of computing / is polynomial -in fixed dimension-). Make the changes 
Zi = Zi tl^" ■ ■ ■ tl!'\ for z G /, in /. Denote by fi each of the rational functions of / after 
the above changes. Next, the Laurent expansion over each rational function, fi, is done in 
the following way: (1) Check if fi contains nondominated solutions computing the Hadamard 
product of fi with /. If fi does not contain nondominated solutions, discard it and set / := I\{i} 
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(termination); (2) if fi encodes nondominated solutions, look for an arbitrary nondominated 
solution (expanding (3) once the first nondominated solution, a, is found, check if there 
exist more nondominated solutions encoded in the same rational function computing /*(/« — 
tj^" ■ ■ ■ t^*"). If there are more solutions encoded in fi, look for them in fi — z"" tj^" ■ ■ ■ t^'"". 
Repeat this process until no new nondominated solutions can be found in fi. 
The process above describes the pseudocode written in Algorithm |2} 

Algorithm 2: A polynomial delay algorithm for solving MOILP 
input : A G Z™^", 6 G Z", C g Z'^^" 

output: The entire set of nondominated [Xe) and efficient {Ye) solutions for MOILPAfi{h) 

Set Xe = {} and Ye = {}. 

Step 1: Compute, f{z), the short sum of rational functions encoding the set of nondominated 

solutions of MOILPA,c{b)- The number of rational functions is indexed by /. 

Make the substitution Zi = Zi tl^' ■ ■ ■ tl^" in f{z). Denote by fi, i E I, each one of the addends in 

Step 2: For each i G /, check fi * f . If the set of lattice points encoded by this rational function 
is empty, do / ^ / \ {i}. 
Step 3: 

while / 7^ do 
for i G / do 

Look for the first nondominated solution, a, that appears in the Laurent expansion of fi. 

Set Xe^Xe^ {a] and Ye^YeU {Ca}. 
Set f,^fi-z'^tr---tr 

and check if / * /j encodes lattice points. If it does not encode lattice points, discard fi 
{I ^ I \ {i}) since fi does not encode any other nondominated point, otherwise repeat, 
end 
end 



Theorem 4.2. Assume n is a constant. Algorithm\^ provides a polynomial delay procedure to 
obtain the entire set of nondominated solutions of MOILPA,c{b) ■ 

Proof. Let / be the rational function that encodes the nondominated solutions of MOILPA,c{b)- 
Theorem ensured that / is a sum of short rational functions that can be computed in poly- 
nomial time. 

Algorithm [2] digs separately on each one of the rational functions fi, i E /,that define /. (Recall 
that/ = ^/,). 

Fix i E P First, the algorithm checks whether fi encodes some nondominated solutions. This 



test is doable in polynomial time by Theorem |1.2[ If the answer is positive, an arbitrary 
nondominated solution is found among those encoded in /«. This is done using digging and the 

n 

Intersection Lemma. Specifically, the algorithm expands fi on the hyperbox ]^[?Tiij, Mjj] fl 
and checks whether each term is nondominated. The expansion is polynomial, for fixed n, since 



the number of terms is polynomially bounded by Lemma 4.1 The test is performed using the 
Hadamard product of each term with /. 
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The process is clearly a polynomial delay algorithm. We use digging separately on each rational 
function fi that encodes nondominated points. Thus, the time necessary to find a new non- 
dominated solution from the last one is bounded by the application of digging on a particular 
fi which, as argued above, is polynomially bounded. □ 



Instead of the above algorithm one can use a binary search procedure to solve multiobjective 
problem using short generating functions. In the worst case, digging algorithm may need to 
expand every nonnegative term to obtain the set of nondominated solutions. Therefore, as 



it is stated in Theorem 4J^, the number of steps to solve the problem can be polynomially 
bounded on M. With a binary search approach, the number of steps to obtain consecutive 
solutions of our problem decreases to a number polynomially bounded on log{M). A binary 
search approach was already used in [TU]. Here, the novelty is that our analysis does not require 
to fix the dimension of the objective space whereas in [TU] it was required. 
The process is as follows. Let M be defined as above. By construction P C [0, M]". We proceed 
by dividing the hypercube [0, M]" into 2" hypercubes of smaller dimensions, and recursively 
repeating the division process over those hypercubes containing at least one nondominated 
solution (until only one solution is included in each element of the partition), whereas those 
hypercubes that at a given stage of the process do not contain nondominated solutions are 
discarded for any further consideration. 

The division process is done by bisecting each dimension. Testing for nondominated solutions 
on a given hypercube (at any stage of the process) is always done using the same tool based 
on Theorem 3.1 That result allows us to construct, in polynomial time in fixed dimension, 
the function h{x) that encodes all nondominated solutions. Moreover, it is easy to see that the 

n 

short rational function that encodes the integer points in the hypercube H = Mi], with 

i=l 

rui, Mi E Q, i = I, ... ,n, is: 



rn[x) 



n rrii Mi 



Thus, the Hadamard product, h{x) * r-H(x) encodes the subset of nondominated solutions that 
lie in 7i; and hence by Barvinok's theory we can also count, in polynomial time, the number 
of integer points encoded by h{x) * r-j^lx) (Lemma 3.4 in |5]). 

The elements in our search space (hypercubes) are organized on a search tree and we use a depth 
first search strategy. Each node is a hypercube containing nondominated solutions. Descendants 
of a given node are hypercubes obtained bisecting the edges on the previous one (parent). It is 
clear that the maximum depth of the tree is 0{logM). The above construction ensures that, 
provided that the set of nondominated solutions is nonempty, finding a first nondominated 
solution can be done testing at most OiT^logM) nodes in the search tree. Since testing a 
node is polynomial, in fixed dimension, this operation is polynomial. Moreover, finding a new 
nondominated solution from a given one is also polynomial. Indeed, it consists of backtracking 
at most 0{logM) nodes until we find a branch containing nondominated points and then we 
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have to explore, at most, 0{2^logM) nodes; or detecting that none of the branches contain 
solutions. 

An illustrative example of this procedure is shown in Figure [T] where can be seen how the 
initial hypercube, [0,4] x [0,4], is divided successively in sub-hypercubes, until an isolated 
nondominated solution is located in one of them. 




Figure 1. Search tree for the problem v — max{{x,y) : x + y < 5,x — 2y < 
2,x + y > 2,x > l,y < 3,x,y & Z_|_} 

The finiteness of this procedure is assured since the number of times that the hypercube [0, M] 
can be divided in 2" sub-hypercubes is bounded by log{M). 
The pseudocode for this procedure is shown in Algorithm [3] 

Algorithm 3: 

Initialization: M = [0, M]" C P. 

Step 1: Let Aii, . . . ,Ai2^ be the hypercubes obtained dividing Ai by its central point. 

i = 1 

Step 2: repeat 

Count nj^4., the number of integer points encoded in rj^.{x) * h{x). This is the number of 
nondominated solutions in the hypercube A^j. 

if TiMi = then 

if i < 2" then i ^ i + 1 

else Go to Step 1 with Ai the next element in the search tree, using depth first search, 
else if n_M. = 1 ( and P fl Mi = {x*} ) then 
I ND = NDU {x*} and i ^ i + 
else 

I Go to Step 1 with M = Mi. 
end 
until ? < 2*" ; 



Theorem 4.3. Assume n is a constant. Algorithm provides a polynomial delay (polyno- 
mially hounded on log{M) ) procedure to obtain the entire set of nondominated solutions of 
MOILPA,c{b). 

Remark 4.1. The application of the above algorithm to the single criterion case provides an 
alternative proof of polynomiality for the problem of finding an optimal solution of integer linear 
problems, in fixed dimension. 
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Assume that the number of objectives, k, is 1, and that there exists a unique optimal value for 
the problem. Applying Theorem \3.1\ ensures that the optimal solution of the problem is found 
in polynomial time, if the dimension n is fixed. 

Remark 4.2 (Optimization over the set of nondominated solutions). In practice, a decision 
maker expects to be helped by the solutions of the multiobjective problem. In many cases, the set 
of nondominated solutions is too large to make easily the decision, so it is necessary to optimize 
(using a new criterion) over the set of nondominated solutions. 

With our approach, we are able to compute, in polynomial time for fixed dimension, a "short sum 
of rational functions" -representation, F{z), of the set of nondominated solutions of MOILPA,c{b) ■ 
This representation allows us to re-optimize with a linear objective, v, based in the algorithms 
for solving single- objective integer programming problems using Barvinok's functions (see e.g. 



35] ) or the algorithm proposed in Remark 4-1 The above discussion proves that solving the 



problem of optimizing a linear function over the efficient region of a multiobjective problem 
MOILPA^cip) is doable in polynomial time, for fixed dimension. 

5. Computational Experiments 

For illustrative propposes, a series of computational experiments have been performed in order 
to evaluate the behavior of a simple implementation of the digging algorithm (Algorithm [T]) . 
Computations of short rational functions have been done with Latte vl.2 [H] and Algorithm 
[T] has been coded in MAPLE 10 and executed in a PC with an Intel Pentium 4 processor at 
2.66Gz and 1 GB of RAM. The implementation has been done in a symbolic programming 
language, available upon request, in order to make the access easy for the interested readers. 
The performance of the algorithm was tested on randomly generated instances for biobjective 
(two objectives) knapsack problems. Problems from 4 to 8 variables were considered, and for 
each group, the coefficients of the constraint were randomly generated in [0, 20]. The coefficients 
of the two objective matrices range in [0, 20] and the coefficients of the right hand side were 
randomized in [20,50]. Thus, the problems solved are in the form: 

(5) max(ci,C2)x s.t. aiXi + ■ ■ ■ + a„ x„ < 6, G Z+ 

The computational tests have been done on this way for each number of variables: (1) Generate 
5 constraint vectors and right hand sides and compute the shorts rational functions for each of 
them; (2) Generate a random biobjective matrix and run digging algorithm for them to obtain 
the set of nondominated solutions. 

Table [T] contains a summary of the average results obtained for the considered knapsack multi- 
objective problems. The second and third columns show the average CPU times for each stage 
in the Algorithm: srf is the CPU time for computing the short rational function expression for 
the polytope with LattE and mo-digging the CPU time for running the multiobjective digging 
algorithm for the problem. The total average CPU times are summarized in the total column. 
Columns latpoints and nosrf represent the number of lattice points in the polytope and the 
number of short rational functions, respectively. The average number of efficient solutions that 
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appear for the problem is presented under eff ic. The problems have been named as knapN 
where N is the number of variables of the biobjective knapsack problem. 



problem 


srf 


latpoints 


nosrf 


mo-digging 


ef f ic 


total 


knap4 


0.018 


12.25 


25.75 


4.863 


4.5 


4.881 


knap 5 


0.038 


31 


62.5 


487.640 


9.25 


487.678 


kiiap6 


0.098 


217.666 


124.25 


2364.391 


7.666 


2364.489 


kiiap7 


0.216 


325 


203 


2869.268 


20 


2869.484 


knapS 


0.412 


3478 


342 


10245.533 


46 


10245.933 



Table 1. Summary of computational experiments for knapsack problems 

As can be seen in Table [T| the computation times are clearly divided into two steps (srf and 
mo-digging), being the most expensive the application of the digging algorithm (Algorithm 
Q. In all cases more than 99% of the total time is spent expanding the short rational function 
using "digging algorithm" . 

The CPU times and sizes in the two steps are highly sensitive to the number of variables. It 
is clear that one cannot expect fast algorithm for solving MOILP, since all these problems are 
NP-hard and #P-hard. Nevertheless, this approach gives exact tools for solving any MOILP 
problem, independently of the combinatorial nature of the problem. 

Finally, from our computational experiments, we have detected that an easy, promising heuristic 
algorithm could be obtained truncating the expansion at each rational function. That algorithm 
would accelerate the computational times at the price of obtaining only heuristics nondominated 
points. 
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